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Abstract. We study a minimal model of traffic flows in complex networks, simple 
enough to get analytical results, but with a very rich phenomenology, presenting 
continuous, discontinuous as well as hybrid phase transitions between a free-flow phase 
and a congested phase, critical points and different scaling behaviors in the system size. 
It consists of random walkers on a queueing network with one-range repulsion, where 
particles can be destroyed only if they can move. We focus on the dependence on 
the topology as well as on the level of traffic control. We are able to obtain transition 
curves and phase diagrams at analytical level for the ensemble of uncorrelated networks 
and numerically for single instances. We find that traffic control improves global 
performance, enlarging the free-flow region in parameter space only in heterogeneous 
networks. Traffic control introduces non-linear effects and, beyond a critical strength, 
may trigger the appearance of a congested phase in a discontinuous manner. The model 
also reproduces the cross-over in the scaling of traffic fluctuations empirically observed 
in the Internet, and moreover, a conserved version can reproduce qualitatively some 
stylized facts of traffic in transportation networks. 



1. Introduction 

The aim of applying methods of statistical physics to the complex behavior of traffic in 
large networked infrastructures is to identify the most important statistical regularities 
and explain the origin of the collective phenomena observed in real systems, such as 
the Internet. As many other complex systems, infrastructure networks can be described 
and studied at different levels of detail, therefore one of the main issues consists in 
recognizing which ingredients are relevant at a given scale, neglecting all redundant 
information. 

In this perspective, many collective phenomena occurring on complex networks can 
be seen as emerging from simple microscopic rules, such as the dynamics of diffusing and 
interacting particles on graphs [1] . Examples of real phenomena that can be explained 
using simple dynamical processes range from avalanches of failures in power-grids [2] , to 
credit contagion in networks of firms [3] and the diffusion of e-mail viruses and spams 
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[5]. The complex phenomena related to traffic, both in transportation [6, 7]. and 
communication networks [8, 9, 10, 11, 12] have been similarly tackled using simplified 
descriptions aimed to characterize their main statistical properties using concepts and 
methods that are typical of non-equilibrium statistical mechanics. The approach of 
statistical mechanics allows to disentangle the complex collective behavior of these 
systems, providing tools to prevent failures and to improve their performances. 

However, the fact of dealing with simple interacting particles systems is sometimes 
perceived as a limit to the richness and variety of the observed phenomena. An 
emblematic counterexample is provided by the discovery that particles condensation 
can occur in one of the simplest classes of non-equilibrium processes on graphs, the 
Zero- Range Processes [13]. The possibility to study these processes analytically, instead 
of resorting to numerical simulations, is of great help to elucidate the properties of 
condensation phenomena in more realistic processes of mass-transfer and traffic on 
networks. 

In the same spirit we propose in this work a simple model of traffic on network 
in the form of a system of particles that are free to hop randomly between nodes but 
are subject to the constraint of forming queues at the nodes. This model was recently 
proposed by the authors as a paradigm to study congestion phenomena on networks 
[14]. Hence, our focus is here congestion in general, as a collective phenomenon arising 
in particles systems in particular dynamical regimes. The mechanisms responsible of 
the phase transition from a free-fiow to a congested stationary state are discussed 
in detail, explaining under which conditions the observed transition is continuous or 
discontinuous. Even though the dynamical process is presented in the very general 
and idealized framework of interacting particles systems, our results have an immediate 
application in the study of several traffic systems, from the Internet to road networks. 
For instance, in the light of our analysis the validity of traffic control strategies in order 
to prevent congestion is questioned. 

The paper is organized in the following way. In Section 2 we present our model 
of dynamical process on graph, discussing its relation with other classes of well-known 
interacting particles systems. A brief account of the phenomenology of the model is 
given in Section 2.2. Section 3 is instead devoted to study the stationary state behavior 
by means of analytical approaches. We provide an iterative method to characterize the 
stationary state on every given graph (Section 3.1) and a mean-field analysis at the level 
of random graphs ensembles (Section 3.2). Finally, in Section 3.3, we discuss the relation 
with condensation phenomena observed in other particles systems. The applications to 
Internet and vehicular traffic are discussed in Section 4. We conclude and indicate some 
possible future development in Section 5. Appendix A and Appendix B are devoted to 
discuss under which conditions a product-measure stationary probability distribution 
exists. Appendix C discusses the case of small queueing capacity, that does not change 
considerably from the large capacity limit considered in the main text. 
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2. The model 

2.1. Definition of the model 

Let us consider a network of N nodes and let v{i) be the set of neighbors of node i. We 
describe particles dynamics as a continuous time stochastic process, in which particles 
arc generated at each node i with a rate pi. Each node is endowed with a first-in first- 
out queue, in which arriving particles are stored. Let Ui be the number of particles in 
the queue of node i. If rii > 0, we assume the following probabihstic hop rule: the 
topmost particle leaves the node at a rate and jumps in the queue of a randomly 
chosen neighbor j e v{i). With probability r]{nj) the particle is rejected by the arrival 
node and remains on the departure one. Otherwise, particles are either destroyed during 
the hopping, with a probability /ij, or, with probability 1 — /ij, enter the queue on node 
J- 

Our model takes an Eulcrian perspective, which focuses on the statistics of the length 
of the queues {ui}, rather than on the trajectories of particles. This perspective allows 
us to disregard the fate of individual particles, by replacing the processes by which they 
are generated and routed to their destination with probabilistic events. 

The long-time behavior of the system is determined by the relation between creation 
and absorption rates. Indeed, if creation rates are much larger than the absorption rates, 
the queues receive more particles of what they can dispose of and the network is rapidly 
overloaded by particles. Such a dynamics leads to the onset of a congested state, where 
queues grow indefinitely. As the external drive imposed by the creation rates does not 
change in time, the system eventually enters a n on- equilibrium stationary state in which 
queues grow with constant velocity. In the non-equilibrium stationary state, a good 
observable is provided by the rate of growth of the total number of particles in the 
system [8], 



where N'it) — '^ini{t) is the total number of particles in the system at time t, 
p — Z^jPi is the average creation rate and r is the observation time. Note that a 
local order parameter, replacing //{t) by ni{t) and p by pj, can be defined in the same 
way. A node is thus congested if, in the stationary state, its average number of particles 
increases with time {pi = {hi) /pi > 0). 

We have already studied this model in the context of packet-transfer based 
communication networks [14], showing the existence of a phase transition horn a free-flow 
regime to a congested state. As we will see later, the character of the phase transition 
depends on the parameters of the model and on the topology of the underlying network. 

A conserved version of the model, without particles' creation and absorption, can 
be defined as well: J\f particles are initially distributed randomly on the nodes, then 
the system is left evolve towards some stationary state in which the distribution of 
the lengths of the queues does not change anymore. The distribution of particles on 



p — lim 



J\f{t + T) -X{t) 

rNp 



(1) 
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the nodes depends on the density of particles p = M /N , on the local hop rates {rj} 
and on the topology of the underlying network. This system presents a condensation 
phenomenon, in which a finite fraction of particles tends to occupy a single node or few 
classes of nodes. 

The model defined here is reminiscent of several other particles systems that 
have been recently studied by physicists and mathematicians, such as the Zero-Range 
Processes (ZRPs), the Misanthrope Processes (MPs) [13], and queueing networks [15]. 
In ZRPs, the hop rate depends only on the number of particles in the departure site, 
whereas in a MP it depends on both departure and arrival sites. The steady-state 
distributions V{ni, . . . jTIn) of ZRPs admit a factorized form, thus these models are 
exactly solvable on every graph. For their simplicity, ZRPs are used as a theoretical test- 
ground for the study of the statistical properties of non-equilibrium systems. Including 
in the hop rate a dependence on the arrival node, the process becomes a MP, that is 
also exactly solvable with factorized steady-state under very mild conditions (see e.g. 
Ref. [13]). The conserved version of our model belongs to this class of processes, with 
hop rates between connected nodes i and j, Uij{ni,nj) = rj[l — r){nj)]/ki. Hence, the 
probability distribution V{ni, . . . , tiat) defining the state of the system factorizes in the 
steady-state in a product measure over single-site distributions Vi{ni) (see Appendix 
A). 

One would be tempted to extend the factorized form to the non-conserved model 
as well. In the Appendix A, wc show under which conditions on the transition 
rates factorization is exact, while Appendix B reports some results on the general 
non-conserved model with particle rejection, from which it is possible to derive an 
approximated mean-field approach. 

The choice of the functional form for the rejection probabilities r]{nj) is motivated 
by the application of the model to several problems from the Internet's dynamics to 
vehicular traffic (see Section 4). In the following, r]{nj) — fjOlrij — n*), that is node j 
refuses particles with a probability f] when it is already occupied by n > n* particles. An 
advantage of using this expression is that as long as the nodes have less than a threshold 
value n* of particles, the model behaves as a "grand-canonical" ZRP, so factorization 
effectively holds. Hence, we expect that in the uncongested state, the product measure 
V{ni, . . . ,nAr) = Hi^il'^i) gives a very good approximation to describe the stationary 
asymptotic regime of the dynamics. This is also confirmed by the absence of correlations 
between queues of different nodes observed in the uncongested phase (see Appendix B). 

2.2. Phenomenology by simulations 

In the previous section we have announced the existence of a phase transition between 
a regime in which particles can freely move in the network and a congested one. This 
can be easily observed simulating numerically the model on any network and varying 
the value of the parameters. Let us consider for the moment the non-conserved model 
with homogeneous parameters /li — /i, Pi — p and ri — l for all i. The simulations are 
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performed on a random regular graph (of fixed degree K = 4) and on an uncorrelated 
random graph with degree distribution P{k) oc k~'^ with 7 = 3. 



2000 




Figure 1. Number of packets as a function of time for an homogeneous network 
of 1000 nodes, degree K — A, without routing procol {rj — 0) 11 — 0.2, in the free 
{p = 0.05) and congested phase {p = 0.25) 

Figure 1 displays two typical time series of the total number of packets A/'(t) in 
the free-flow and congested phases, obtained varying the creation/ absorption ratio p/ fi. 
In the free-flow phase, the number of particles fluctuates around a stationary value 
that is much smaller than the network's size A^, i.e. most of the nodes are empty. In 
the congested phase, instead, the number of particles waiting in the queues constantly 
increases in time, with the lack of stationarity in the number of particles. In this example 
the congested phase appears around Pc = fi, that seems trivial if we notice that a single- 
queue server with constant arrival rate p and departure rate /i becomes overloaded 
exactly at = 1. In fact, the threshold Pc = fi is correct only on homogeneous 
networks (random regular graphs, regular lattices, etc.), not in the presence of large 
degree fluctuations. 

The difference between homogeneous and heterogeneous networks, as well as the 
role played by particles rejection, becomes evident looking at the behavior of the 
congestion parameter p{p). In Fig. 2 (left) we report p{p) for a random regular graph in 
two different situations: with small rejection (f/ = 0.1, n* = 10) and with strong rejection 
{f] = 0.9, n* = 10). The same plots for a scale-free network with 7 = 3 are reported 
in Fig. 2 (right). Let us focus on the curves obtained with small particles rejection. 
The homogeneous network becomes congested with a continuous phase transition taking 
place at Pc — /i, whereas for the heterogeneous one the congested state appears much 
earlier at p <^ yU. In both cases, the effect of the increasing of f] is that of changing 
the nature of the phase transition from continuous to discontinuous with hysteresis. 
The continuous and discontinuous transitions (coming from lower p) are located at the 
same point pc — ^ on homogeneous networks; on heterogeneous networks, instead, the 
discontinuous critical point is shifted towards larger values of the creation rate (with 
respect to the continuous critical point). 
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Figure 2. Left: Transition curves p{p) for a random regular graph of size N = 10000, 
H = 0.2, fj = 0.1, fj = 0.9. 

Right: Transition curves p{p) for an uncorrelated scale free graph with 7 = 3, k„iin = 2 
oi N = 3000 nodes, /i = 0.2, f] = 0.1, f] = 0.9 for n* = 10. For 77 = 0.9 the system 
shows hysteresis in both the homogeneous and heterogeneous case. 



The theoretical understanding of these congestion phenomena depending on the 
hop rules and the underlying networks, as well as their potential application in the 
study of real traffic problems are the subjects of the next sections. 

3. Analytical approach to the stationary state 

3.1. Iterative equations on single graphs 

In the conserved model, the single graph analysis is simple because factorization is exact, 
and the result on any given graph is exposed in Appendix A. The situation is different 
when particles are created and destroyed. Nevertheless, correlations between nj's on 
different nodes are hardly detectable in numerical simulations, as shown in Appendix 
B. This, and the fact that the interactions are local, makes it reasonable to work within 
an factorized approximation, i.e. V{ni, . . . .n^) — ni^il"-*)' ^^"^ derive some local- 
recurrence equations that can be solved numerically in polynomial time on every graph. 
In this way, one can study very general assignments of parameters and node functions 
{Pi,/ii,ri,r7}. 

Here we present a derivation of these iterative equations that is based on a simple 
detailed balance argument for the single-node dynamics in the factorized approximation. 
A more precise derivation is presented in Appendix B, where the same equations are 
obtained by means of a series of approximation starting from very general exact results. 

In the single-node description, the transition rates for the queue length of node i 

are 

w{n, ^ + 1) = + (1 - - r/(n,)) "^^'^^ ~ (2) 
where v{i) is the neighborhood of i. A reasonable choice for the rejection probability 
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is rj{ni) = f}9{ni — n*), where node i is congested as soon as > n*. Imposing the 
detailed balance, the distribution Vi{ni) turns out to be a combination of exponentials, 
depending on rii > n* or not (see also Appendix B). 

We expect three different behaviors: free nodes with exponentially decreasing dis- 
tribution, congested nodes with exponentially increasing distribution (not normalizable) , 
and unstable nodes, whose distribution is peaked around n*, that we call fickle nodes. 
Assuming the validity of the double-exponential single-node distributions, it is easy 
to show that the stationary state can be described in terms of only two local quan- 
tities: Qi — Prob{ni — 0} e [0,1], the probability of empty queue at node i, and 
Xi — Prob{xi = 1} e [0,1], the probability that node i has more than n* particles 
{xi = e{ni - n*)). 

A congested node has always Ui > n*, thus Xi = 1 S'-iid qi = 0. Using Eqs. 2, the 
average rate of increase of Ui in the stationary state is given by 

= p. + (1 - - fiX^) E "-"^^ - ^-^^ E [1 - ^XA- (4) 

jev{i) jev{i) 

It is straightforward to verify that > for congested nodes. 

Free nodes always have rij < n*, thus Xi = 0- addition Ui = otherwise they will 
eventually become congested. Imposing this condition on Eq. 4, we get 

P» + (1 - ^i) Eje^(») ""'^"^'^ 
qi = Qiix,q) = ^ T^^TV ^ ■ 

The case of fickle nodes is more tricky because their probabihty to be empty is expected 
to be exponentially small with n*, i.e. oc e~"»*, and vanishes only for n* —>■ oo. For the 
sake of simplicity we consider such a limit, that is approximately correct when n* ^ 1. 
Hence, imposing fii = and neglecting <^ 1, we find an expression for Xi, 



Xi = Ciix, 0) = ^ 
V 



(6) 



(1 fJ'i) ^j€v{i) k 

It is easy to check that for a congested site Ci(x, q) > I and Qi{x,q) < 0, therefore the 
previous expressions can be written in the more compact way 

Xi = max{0,min[l,Ci(x,g)]}, .^x 
qi = max{0, min[l,(5i(x,g)]} ■ 

These self-consistent equations can be solved iteratively on any specific graph. 

If a fixed point of Eqs. 7 exists, the global congestion level can be measured by 
the order parameter p — ■^Yliii'^i) i where (rij) is the average growth rate of queue i 
computed on the fixed point values. 

As an example. Fig. 3 shows the diagram p{p) obtained solving Eqs. 7 (in the limit n* — > 
CO ) on a random regular graph (left) and a scale- free network (right) for an homogeneous 
choice of the parameters {pi = p, pi = p = 0.2, and f] = 0, 0.25, 0.5, 0.75, 0.9, 1). In the 
random regular graph (left), the case f] = 0, corresponding to a ZRP with creation and 
absorption of particles, presents a clear signature of a continuous phase transition from 
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Figure 3. Behavior of tlie congestion parameter p(p/ fi) for a random regular grapli 
of size N = 10'^ and degree K — 4 (left) and a scale-free network of size N — 10'^ 
and exponent 7 = 2.5. The symbols represents the behavior obtained by solving 
numerically the iterarive equations 7 for increasing values of p and rj — (black 
open circles), 0.25 (red full circles), 0.5 (blue squares), 0.75 (green crosses) 0.9 (violet 
triangles), 1 (black crosses). Solving the Eqs. 7 for decreasing values of p, we find 
instead the corresponding dashed curves. 



a free-flow regime to a congested phase. When f/ > 0, we expect hysteresis phenomena, 
thus we first consider the solutions obtained for increasing values of p. More precisely, 
we find a fixed point for a value of p, we slightly change p and we iterate the Eqs. 7 
starting from such solution until we find another fixed point. Increasing p from zero, the 
congested phase appears abruptly, with a discontinuous behavior, at the same critical 
rate Pc = A* (independently of the value of f]). On the contrary, decreasing p from the 
congested region, the systems undergoes a continuum transition to the free-flow phase 
(dashed lines) whose position decreases increasing f/. 

In the scale-free network (right) the critical value Pc ^ for f/ = 0, but it is shifted 
towards higher values as soon as f/ > 0, as already observed in Fig. 2 and found in 
[10, 14]. It is worthy noting that in scale- free networks the position of the transition 
depends on the value of f] and the curve p{p) changes convexity for increasing values of 
f]. Finally, for f/ = 0.75 the transition occurs in two-steps, first a continuous transition 
then a discontinuous one at slightly larger values of p. As expected, decreasing p from 
the congested phase we observe hysteresis (dashed lines). We have verified performing 
the same calculation on other networks that the double transition is not always present 
and depends on both the value of fj and the tail of the degree distribution. 

Though our results on single graphs reproduce the phenomenology of congestion 
observed in previous numerical simulations [10], they also pose many new questions 
about the nature of the phase transitions and the mechanisms behind them. For this 
reason, we have developed a complementary mean-field analysis at the level of random 
graphs ensembles that is able to shed light on all these points. 
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3.2. Mean-field analysis at the ensemble level 

We consider uncorrelated random graphs with degree distribution P{k), so that rik 
represents now the average queue length of nodes in classes of degree k. For the sake 
of simplicity, we will examine the simple homogeneous case Pi — p, /li — /i, ri — 1 and 
77(71) = fj9{n — n*). We define Qk — P{ni = 0\ki = k} as the probability that a node 
of degree k has empty queue, and Xk = vPi'^'i ^ n*\ki = k} as the probability that a 
node of degree k refuses particles (it should be noticed that now x has an 77 factor). The 
mean-field transition rates for nodes with degree k are 

k 

Wk(n ^7i + l)=p+(l- //)(! - q)-{l - fie(n - n*)) 
Wk(n^n-l)^e(n)(l-x), (8) 

where z is the average degree, q = qkP{k) and x = Ylk ^XkP{k). The average queue 
length (uk) follows the rate equation 

(hk) = p + (1 - //)(! - g)-(l - Xk) - (1 - qk){l - X). (9) 

z 

Note that summing over k and dividing by p we obtain a measure of the order parameter 
pip)- 

Since fik depends linearly on k, high degree nodes are more likely to be congested, 
therefore, for every p, there exists a real valued threshold k*{p) such that all nodes 
with k > k* are congested whereas nodes with degree less than k* are not congested. 
Congested nodes {k > k*) have qk = and Xk = V- The probability distribution for the 
number of particles in the queue of free nodes with degree k < k* can be extracted by 
calculating the generating function Gk{s) — Vk{nk — n)s'^ from the detailed balance 
condition Wk{nk + 1 ^ nk)Vk{nk + 1) = w{nk — > 71^ + l)Vk{nk), that we assume to hold 
in this approximation (see also Appendix B). The generating function takes the form 

G>M-«J^^^+ , 'r^"' , I (10) 

1-akS l-{ak-bk)sj 

corresponding to a double exponential, where ak = [p + (1 — A*)f (1 — q)]/[^ — x] aiid 
bk = fj[{l — A*)|(l — ^)]/[l — x]- From the normalization condition Gfc(l) = 1 and the 
condition fik — 0, we get expressions for qk, Xk, 

-1 



qk 



^ ^k 



1- Qk 1- ak + b 



(11) 



p - (1 - qk){l - x) /-i^x 

and, finally, for q, x- 

The value k* is self-consistently determined imposing that nodes with k — k* are 

marginally stationary, i.e. fik* — with qk* — 0, Xk* — V, that translates into the 
equation 

k* = rz. (13) 

(l-/.)(l-77)(l-g) ^ ' 
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The set of closed equations for q, x can be solved numerically for any degree distribution 
P{k) and p{p) can be accordingly computed. 

3.2.1. Homogeneous Networks - The equations for q and x simplifies to a single 
equation when all nodes have the same properties, and in particular the same degree 
{ki = K, \fi). On these networks, the mean-field behavior can be trivially studied 
for any value of n*, but we consider as an illustrative example the limit n* — > oo. 
Only two solutions of the equation relating q and x are possible: the free-flow solution 
(p = 0) with q = 1 — p/ fi and x = that exists for p < fi, and congested-phase 
solution, where all nodes have rii oo, i.e. x = V q = 0. The latter solution has 
p = n/p = 1 — {1 — ff)p/p and exists for p > {1 — fi)fi. There is a simple argument 
explaining the law 1 — {1 — f])n/p, that corresponds to the situation in which the whole 
system is congested. In the infinitesimal interval of time dt, pN particles are created, 
and particles try to hop. The nodes are congested, so a fixed fraction p{l — fj) of 
them is absorbed, and the expression of p{p) follows from its definition (Eq.l). 

The behavior of the congestion parameter with both the continuous and 
discontinuous transitions to the congested state is plotted in Fig. 4 for fj = 0.25, 0.75. 
The corresponding phase diagram, reported in the inset of Fig. 4, shows that in the 
interval p G [(1— f/)/i, p] both a congested- and a free-phase coexist. We find an hysteresis 
cycle, with the system that turns from a free phase into a congested one discontinuously 
as p increases and crosses p = p, and it reverts back to the free phase only at p = (1— f/)/i 
as p decreases. 




Figure 4. Behavior of tlie congestion parameter p{p/ fi) for a random regular network 
obtained theoretically for rj = 0.25, 0.75. Inset: phase diagram for the same graph. 



3.2.2. Heterogeneous Networks - In the case of heterogeneous networks the equations 
for q and x have to be solved numerically. For instance, in Fig. 5 we compare 
the theoretical prediction (full line) for p{p) in a scale-free network with results of 
simulations (points). The agreement is good, the theoretical prediction at the ensemble 
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level confirming the scenario already observed in the simulations of Section 2.2 and 
in the single-graph analysis of Section 3.1. The curves are obtained for fi = 0.2 and 
n* = 10, but the behavior does not qualitatively change for different values of these 
parameters. The dependence on fj brings instead qualitative changes. Increasing f/ from 
0.1 to 0.9, the transition becomes discontinuous and pc increases. 




Figure 5. p{p) for an uncorrelated scale-free graph (-P(fc) cx fc ^, fcmin = 2, 
= 110, N = 3000), n = 0.2, n* = 10 and fy = 0.1 and f] = 0.9, from 
both simulations (points) and theoretical predictions (lines). Hysteresis is observed 
increasing (black curve and points) and then decreasing (red curve and points) p across 
the transition. 



The main difference with respect to homogeneous networks is that on heterogeneous 
networks, not all nodes become congested at the same time. The rate p at which a node 
becomes congested depends on its degree, the hubs being first. The process governing the 
onset of congestion and the effects of the rejection term can be understood in the limit 
n* oo, that simplifies considerably the calculations without modifying the overall 
qualitative behavior for sufficiently large n*. We refer the reader to Appendix C for a 
discussion of the main differences emerging in the case of low values of n*. 
We have to solve in the limit n* —>■ oo the self-consistent equations for x and q. In this 
limit, uncongested nodes have a^. < 1, hence Xk ^ and = 1 — a^. All nodes with 
degree k < kp, where kp = max(fc*(l — f]),kmin), are free from congestion. Congested 
nodes have qk ^ and Xk = V (for k > k*). The fickle nodes are those with kp < k < k* 
and they have Xk = ^ ~ Using this classification, we get a first expression for x, i.e. 



Xi 



k 

k_i 
k 



E i-¥ -p(k) + nE-m- (14) 

k=kf 

Eq. (13) provides a further relation between q, x and k*. We eliminate q using its 
definition which leaves us with another expression for x. 



X2 



1- — |l + Ap-fi + [{1 + Ap - + AABpY^^'j (15) 
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where A = z/[k*{l - fj){l - //)] and B = Y!l=kmiu ^ ~ ^(^)- determine x we 
have to solve the imphcit equation xi — X2- 

In Fig. 6 we plot the difference Ax = Xi ~ X2 vs. k*, for 77 = 0.1 (left) and 0.9 (right) 
and different values of p on a scale-free graph. The zeros of Ax{k*) correspond to the 
only possible values assumed by k*. For small rejection probability {fj = OA in Fig. 6), 
there is only one solution k*{p), which decreases from +00 when increasing p from 0. 
The value Pc at which k*{pc) = kmax is the critical creation rate at which largest degree 
nodes become congested. At larger p, k*{;p) decreases monotonously until eventually 
all nodes are congested when k*{p) — kmin- Hence for low values of 77, the transition 
from free-flow to the congested phase occurs continuously at the value of p for which 

k*{p) = kmax- 

At large fj {fj — 0.9 in Fig. 6), the scenario is more complex. Depending on p, the 
equation can have up to three solutions, kl{p) < /i;|(p) < k^{p). It is easy to check 
that only kl and k^ can be stable solutions. For p <^ 1 there is only one solution at 
^Kp) ^ kmax, corresponding to the free phase. This is thus the stable solution for p 
increasing from zero. As p increases, another solution kl{p) < k^{p) can appear, and 
kl{p) moves towards lower degree values. Three situations may occur: 

i. The solution kl{p) disappears before reaching k^ax- Then kl{p) becomes the stable 
solution, and the congested phase appears abruptly. However, given the shape of 
the function Ax (see Fig. 6), when this happens kl{p) — > and in particular we 
expect klip) < kmin, so that above the transition the whole network is congested 
and follows the law p{p) = 1 — (1 — 77)^. 

ii. The solution k^ip) crosses kmax and exists until it reaches kmin- Then the congested 
phase emerges continuously and the network is only partially congested (i.e. only 
the nodes with k > kl{p)). The order parameter grows until it reaches the curve of 
complete congestion p{p) = 1 - (1 - 77)^ (A;|(p) < /oWn)- 

Hi. The solution k^{p) crosses kmax but disappears before reaching kmin, and kl becomes 
the stable solution. In this case the congested phase appears continuously (only 
high-degree nodes are congested), but at some point another transition occurs that 
brings the system abruptly into the completely congested state. 

The possible situations are summarized in Fig. 7, where we have sketched the 
corresponding behaviors of the order parameter p{p). 

The scenario at points i.-ii. is exactly that of Fig. 5, while a signature of the 
double-transition can be observed in Fig. 3 (right). In general, the exact phenomenology 
observed in numerics and simulations depends strongly on the tail of the degree 
distribution, i.e. on the graph ensemble considered. 

Note that in case of discontinuous transitions, the presence of an hysteresis 
phenomenon is associated to the stability of the two solutions kl{p) and kl{p). For 
instance, in case ii or Hi, we start from the free-phase at low p, the system selects the 
solution k^{p) and follows it upon increasing p until the solution k^{p) disappears. On 




Figure 6. The zeros of Ax(p) vs. k* define the threshold degree for the onset of 
congestion in a network. The picture refers to a scale-free random network with 
7 — 3.0, kmin ~ 2 and N ~ 3000 {k,nax = HO), and different values for fj — 0.1 
(left) and 0.9 (right) and p. The solution k\{p) in the right panel falls outside the plot. 



the contrary, starting from the congested phase (large p) the system selects the solution 
klij)) and remains congested until this solution disappears (see inset of Fig. 5). 

In Fig. 8 we can see the solution k*{p) for the same graph of Fig.5, with f/ = 0.7: 
at pi, when k* = kmax, the system becomes congested in a continuous way, at there 
is a discontinuous jump to higher values of congestion, while above p4 the network is 
fully congested and finally, coming back to p2 there is a jump to a less congested state. 
Between p2 and p^ there is coexistence of high and low congested states with hysteresis. 




Figure 7. Increasing fj, the congestion parameter p{p) develops a discontinuous 
transition. Here we report the case of the graph of Fig. 6. For rj — 0.75, we have 
first a continuous, then a discontinuous transition. 



In summary, the system can show a sort of hybrid transition: a continuous transition 
to a partially congested state followed by a discontinuous one to a (almost) completely 
congested one (see Fig. 7). 

3.2.3. The General Phase-Diagram On heterogeneous random graphs, the behavior of 
the system in the plane (f/,p) depends in a complex way on its topological properties, 
such as the degree cut-off and the shape of the degree distribution. For this reason 
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Figure 8. The solution k*{p) for the scale-free graph of Fig. 5, with f] = 0.7. At pi, 

k* — kmax, and the system becomes partially congested in a continuous way. Between 
P2 and P3 there are three solutions, two of them are stable. Increasing p, the system 
jumps suddenly to a more congested state at ps, whereas decreasing p, the system 
jumps to a less congested state at p2- Above pi the system is completely congested. 



the precise location of the critical lines, separating different phases, can be determined 
only numerically using the methods exposed in the previous section. In the following, 
we give a qualitative description of the general structure of the phase diagram in the 
limit n* — > oo, then we substantiate the analysis reporting an example of phase diagram 
obtained numerically for the same networks ensemble of Fig. 5. 

A first important region of the space of parameters is the one in which a completely 
free solution exists, i.e. k^ax < kp. This solution is characterized by ^ = 1 — p///, x = 
and p = 0. Prom the expression for = computed in kmax we find that this happens 
as long as p < Pco with 

Pco = , ^ • (16) 

Note that this region does not depend on the rejection probability f^, because rejection 
affects only congested nodes. 

The transition takes place when the maximum degree nodes first become congested, 
i.e. k* = kmax- Since hk* = 0, qk^^^ = and Xkma. = V, we get fi:om Eq. 9 a first 
expression for Pc = 1 — X — hm^(^\ _ _ _ ^) Now computing p averaging Eq. 9 
and imposing p = 0, we find a second expression for — //(I — q){l — x). Eliminating 
q from these two equations, we find the critical line 

where % = Y.k>kF (l " • Below this fine (dotted line in Fig. 9) the 

system is not congested {p — 0), even if in the region Pco < P < Pc{v) higher-degree 
nodes are unstable {kp < kmax ^ k*). 

It is possible to show that Pciv) attains its maximum in fjc = 1 — jr^ where Pcmax — 
nhmin^ where kp — kmin a-^d so above this point the curve is constant Pc{fj > fjc) — PciVc)- 
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Figure 9. {r],p) phase diagram for the uncorrelated scale- free graph of Fig. 5. 



The transition line pdv) corresponds to the point pi in Fig. 8, calculated for all 
values of fj. We can calculate the two curves P2{v), Psiv) as well, in order to get the 
points at which there are discontinuous jumps in the congestion parameter p(p). 

Looking at Fig. 9 we can distinguish three points A, B, C dividing the phase 
diagram into different regions: 

i. Below fjA we have a continuous transition to a congested state increasing p above 



ii. Between fj^ and fjB the transition is continuous at pi. Then, increasing p above 
P3, there is a discontinuous jump to a more congested state. Coming back to lower 
values of p, there is a discontinuos jump to a less but still congested state at p2, 
and the system eventually becomes free below pi in a continuous way. 

Hi. Increasing p in the region between fis and fjc, there is a smooth transition from 
free-flow to a congested state at pi, and a sudden jump to a more congested phase 
at ps', but, this time, by decreasing p from the congested state, the transition to 
the free phase is discontinuous and located in p2. 

iv. For fj > fjc the transition is a purely discontinuous one with transition points p2 



Increasing p above the transition, at some point Pciiv) the system becomes 
completely congested. For p > Pciiv), the order parameter follows the curve p = 
1 — /Li(l — ri)/p. This happens for p > Pdiv) = (1 — ~ ~ f^)kmin/z), where 
k* < kmin, q = 0, x = V- 

These calculations show that the phase diagram crucially depends on the tail of the 
degree distribution. In scale-free networks k^ax scales with the network's size N as 
with u; — 2 (structural cut-off) or a; = 7 — 1 (natural cut-off). Accordingly the critical 
line depends on the system's size, pc oc N~'^ . The only region that does not depend on 
kmax is that after the maximum of the curve {fj > fjc)- 



Pi- 



and P3. 
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3.3. Relation to condensation phenomena 

We have discussed the stationary dynamics of the non-conserved model, but it is still 
not clear if the mechanism triggering the congestion phase transition is the same causing 
condensation in ZRP [13]. 

The steady-state properties of the conserved model (p = 0, // = 0) are exactly 
solvable once we have factorized the distribution on the nodes. Imposing the detailed 
balance we get (the calculation with rj 7^ 1 is in Appendix A), 

^" \ P.(0)[^]"(l-77)— * n>n* ^^^^ 

where Vi{0) is determined from the normalization condition on node i and A is a constant 
independent of the node that is fixed by the total number of particles J2i '^i — Their 

explicit expression is not necessary for the purposes of the present analysis. We can 
interpret s = 1/A as a chemical potential that fixes the number of particles in the 
network, and study the "grand-canonical" generating function 

N 



n 



1-ski ' 1 - (1 - fj)ski\ ■ ^"^^-^ 

In the grand-canonical formulation, the particle density can be written asu = jj^^^^^^^ 
[13]. We focus now on the limit n* — > oo, for which the particle density becomes 
i/(s) = i The sum is defined for values of s smaller than the smallest pole of 

the argument, i.e. s < 1/kmax, i-e. for A > kmax- So as long as A > kmax the integral 
converges and the system is able to allocate the corresponding density of particles on 
the network. This is not possible when A — > kmaxi because a pole appears at kmax- 
We follow a recent approach by Noh [16] and isolate the contribution of the maximum 
degree nodes vm{A) = ;^ Y^f^^^- Grouping together nodes of the same degree, we find 
in the continuous approximation 

/kmax P(k)k 
dk^^ (20) 

where P{k) oc k~^ is a power-law degree distribution. In random networks the maximum 
degree scales with the system size as k^ax ~ f^maxN^^^, where to depends on the 
generating algorithm (e.g. a; = 2, 7 — 1), it is thus natural to rescale the variables 
in Eq.20, A aN^^'^ and to use x = k/A, obtaining 

kmax/A X^~^ 



v[a) — vuio) + I dx 



I A 1 -a^ 



= VM{a) + O (tV-^) - O (ivV) log(l - Kmaxia) (21) 
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Figure 10. Plot of the function v{p) in Eq. 23. The divergence occurs at 
p = Pc = f^/ll-i' + (1 — l-i^kmax/z]. The inset shows that the divergence is triggered 
by the behavior of kmax but nodes of increasing degree contribute to increasingly large 
portions of the global density. 



where we have estimated the divergence of the integral in the two extremes of integration 
for Kmax/d — > 1 and Kmin/d 0. Neglecting the second vanishing term and replacing 
the ratio Kmax/O' — NumH^ + Ni/m), we get 

u{a) ^ UM{a) + 0{N^) log(l + NuM{a)). (22) 

The last term in the r.h.s. vanishes in the infinite size limit, whereas the first term is 
finite. This means that for large systems condensation takes place and a finite fraction 
of the density is stored in the maximum degree node. Decreasing further A < kmax, the 
condensation occurs on all nodes of degree k > A, and the condensate is not localized 
on a given node. This is an important difference between condensation induced by 
topological heterogeneity [16] and standard condensation in homogeneous systems [13]. 

The extension to the original non-conserved model is not straightforward, for the 
absence of a factorized solution, but a qualitative understanding can be obtained using 
a mean-field approximation. In this way, we can compute the particle density from 
the generating function (Eq. 10) following the approach of Section 3.2. In the simple 
case n* — > oo the sum can be evaluated analytically, but in this limit the sum diverges 
as soon as the first fickle nodes appear (because their queues have distribution peaked 
around n* that is taken to infinity). The curve exists when kp < kmax, so x = and 
1 — q — p/ /I, thus 

^ p+(l-/x)^^ 



Z IJ, 



u{p) diverges continuously (see Fig. 10). The 



At the critical point Pc — —m — ^ / : 
divergence is triggered by nodes of maximum degree, with the same mechanism of 
condensation, and increasing p above the transition the set of nodes on which particles 
accumulate grows including nodes of lower and lower degree. In Fig. 11 we report the 
behavior of the density of particles as a function of p for systems with finite n*. The 
data, obtained numerically simulating the system for fj = 0.9, are compared with the 
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Figure 11. Plot of the density profiles ^{p) with {rj = 1 and several n*) and without 
routing protocol from numerical simulations run on the same scale-free graph as in 
Fig. 6. 



theoretical prediction, obtained solving the man-field equations for finite n* = 2, 5. The 
agreement is reasonably good and as expected the density increases for larger values of 
n* because more particles can be stored into the queues. 

4. Applications to real traffic systems 

We have put forward a simple model of non-interacting particles moving randomly 
and forming queues on networks, and used it to understand the collective behavior that 
induces congestion phenomena on networks. A possible application is to describe packet- 
transport in information networks, or vehicular traffic on transportation networks. In 
the following we briefly discuss these two problems and the properties that can be 
correctly predicted by simple models inspired to the present one. 

4.I. The Internet at routers level 

Congestion in packet-based communication networks means that the network, or a part 
of it, is not able to process all arriving information, nodes become overloaded, and 
this implies a global slowing-down of the system. Congestion phenomena have been 
observed in wireless networks [17], in multimedia networks [18], and, more importantly, 
in the Internet [19]. The first identified Internet's congestion collapse dates back to 
October 1986, when data throughput from LBL to UC in Berkeley suddenly dropped 
from 32 Kbps to 40 bps. After that initial event, traffic congestion continued to threaten 
Internet's practitioners, because of the impossibility of constantly monitor and supervise 
large portions of the Internet, and clearly identify precursors of a congestion event. For 
these reasons, understanding congestion phenomena in packet-based communication 
networks has become a subject of intense interdisciplinary research [4], with many 
contributions from statistical physicist community, particularly after the works by 
Takayasu and collaborators [9], in which the evidence of a phase transition from a 
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free-flow regime to a congested phase depending on the load level was reported. 

In information networks, like the Internet, packets of information arc created at 
some nodes and then forwarded nodc-by-node until they reach their destination. The 
packets are dispatched by a routing protocol that tries to minimize the traveling time, 
taking into account information about the distance and the local traffic. Most of the 
theoretical models for such a dynamics are based on a shortest-path routing protocol, in 
which packets follow the shortest-path between a given pair of nodes (where packets are 
created and destroyed). Echenique et al. [10] have found in numerical simulations 
that the nature of the transition depends on the type of routing rules: in case of 
purely topological routing (e.g. along the shortest paths), the congested phase appears 
continuously, whereas the transition is discontinuous if some traffic-aware scheme is 
considered (e.g. delivery packets preferentially to uncongested nodes). Other works have 
considered more general forms of routing rules, with particular emphasis on optimization 
strategies to improve network performances [11]. 

The microscopic dynamics of our model seems quite different from these realistic 
processes, but it contains all features that are relevant to characterize the free- 
flow/congestion transition in information networks, reproducing the collective behaviors 
already observed in the literature. The fact that the absorption of packets occurs only 
when packets move, not when they are waiting in the queues, mimics the behavior 
of real packets of information that leave the network when they reach a destination 
node. On the other hand, random walks and shortest-path routing have very different 
statistical properties. The visiting probability of a node i in the shortest-path routing 
is proportional to the betweenness centrality of that node, and thus scales non-linearly 
with its degree, whereas in the random walk protocol the relation is linear. In order to 
accommodate for this statistical feature, and make the routing process more realistic, 
one could consider degree-biased random walks [21]. Another important ingredient 
of the model is the presence of a rejection probability r]{ni), that reproduces the 
congestion avoidance scheme elaborated by computer scientists for the Internet [22]. 
This class of algorithms are based on a feedback mechanism that relies on the exchange 
between routers of Acknowledgement signals (ACKs) carrying information on the local 
level of traffic. When the round-trip-time of ACKs sent in a given direction becomes 
too large, the node decreases the rate with which packets are forwarded in such 
direction. Therefore, like in the present model, in the Internet congested nodes have 
a lower probability to receive packets. Such a scheme is useful to retard the onset 
of the congested phase, but it introduces a cooperative behavior that can lead to a 
discontinuous transition. 

In order to characterize the kind of congestion phenomena that could be observed 
on more realistic topologies than random graphs, we have studied our traffic model 
on an Internet's map at the routers level obtained by the CAIDA group of Internet's 
measurements [23]. The map counts N = 192244 nodes, maximum degree k„iax = 1071 
and average degree z — 6.3. It has a degree distribution that is well-fitted by power-law 
(P(/c) oc k~'^ with 2 < 7 < 3). Figure 12 (left) shows the behavior of the congestion 
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Figure 12. Transition curves p{p) for the CAIDA network of routers obtained by 
siniulations(points) and theoretical predictions (lines) with fj = and fj = 1. We have 
also reported the corresponding fraction of congested nodes (dotted lines) 

parameter p{p) on the routers network obtained running simulations of our traffic model 
with and without congestion-aware routing protocol. In both cases the congested phase 
emerges continuously, but with very different behaviors. At low p values the protocol is 
able to reduce the congestion level, but at larger values of p it is no more effective and 
the level of congestion starts increasing much faster for fj = 1 than for f/ = 0. These 
results are confirmed by the numerical solution of iterative equations Eqs. 7 shown in 
Fig. 12. In Fig. 12, we report also the behavior of the fraction of congested nodes 
(dotted lines) for both f/ = and 1. Comparing the two sets of curves we observe 
that in the absence of the traffic-aware routing protocol the congested phase is mostly 
concentrated on a subset of nodes (i.e. the most highly connected ones), whereas for 
fj > the congestion is homogeneously spread over all congested nodes independently 
of their topological properties. 

Another characteristic trait of the Internet is the scaling crossover empirically 
observed in the single-node's traffic statistics. More precisely it consists in measuring 
the scaling of traffic's fluctuations cr^(n) with respects to the average traffic on a node 
(n). In Ref. [25], it was found that traffic on nodes at the level of the Abilene Internet's 
sub-network presents two distinct regimes. At low traffic level the fluctuations are 
linear in the average traffic cr^(n) oc (n), whereas at higher levels of traffic they found 

oc (n)^ (with prefactors possibly depending on the node's degree). 
In the present case, for a free node i of degree k, the fluctuations of the queue's length are 
easily obtained within the mean-field approximation (with n* oo). The generating 
function gives Gk{s) = 1/(1 — a^s) with ak < 1, and using G'j.{l) = (uk), we get 
al = (uk) (1 + (rife)). At low traffic levels (uk) <^ 1, thus al ~ (rife). Instead when the 
traffic increases, (rife) becomes larger than 1 and cr^ ^ (nk)"^. 

We have computed the scaling behavior on the realistic CAIDA's map, using 
numerical simulations of our model. As expected, the model correctly reproduces this 
important statistical feature, as evidenced in Fig. 13, in which we display the fluctuation 
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Figure 13. Scaling of tiie fluctuations respect to the mean from simulations onto the 
CAIDA network. 

scaling obtained averaging over the whole system (i.e. average over the degrees). 
4-2. Traffic in transportation networks 

Congestion phenomena are observed in transportation networks as well [24, 26, 27], 
even if the dynamical features can be very different from those defined on information 
systems. A first difference is that vehicles flow along roads that are represented by links, 
whereas the nodes are just the intersections between them. This is usually accomodated 
by means of a dual approach, in which nodes and links are exchanged. More importantly, 
road networks have a very peculiar structure with significant constraints imposed by the 
real space embedding that influence both the large-scale topology (e.g. minimization 
of the distance, planar structures), and the single- node properties (road have finite 
capacity). All these features are expected to have an effect on traffic and thus on the 
appearance of congestion. 

Traffic in urban road networks is usually studied by means of agent-based or 
fiuidynamical models, with a large quantity of realistic ingredients [27]. Nevertheless, 
at a purely statistical level, a simple model of random walkers with queues can already 
provide a correct qualitative description of the collective behavior of the system. 

We represent schematically vehicular traffic on a road network considering random 
walkers on a planar graph. The fact that vehicles proceed along roads in an ordered 
way is mapped into a queueing process on the nodes of the dual network. For simplicity 
we fix the number of particles (or the density u = JV/N), ignoring the agents that 
may enter or leave the system. Since roads have a finite capacity, we fix the maximum 
number n* of particles that a node can store in its queue; and in order to avoid local 
overloads, we assume complete rejection when the queue has reached its capacity, i.e. 
ri{ni) = 1 for > n*. Moreover, when the density of cars on a segment of a road is 
not too high, cars use to travel at a constant speed, that can be mimicked by letting 
the particles move as non- interacting random walks until a density threshold ui. At 
high density, the flux becomes constant because cars form queues. Thus we assume that 
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Figure 14. Outgoing flux u from a road (noad) as a function of the number of 
particles n. Until ni it is proportional to n, and then is constant. The maximum 
number (capacity) is n* . If the destination node has n* particles the flux is zero. 



above ni, particles form queues as well with constant outgoing rate. In summary, the 
flux u out of a road (node) is (see Fig. 14): 

i. u{n) = n if n < rii and the destination node has less than n* particles. 
a. u{n) = ni ii n > ui and the destination node has less than n* particles. 
Hi. u{n) = if the destination node has n* particles. 

Vehicular traffic is usually studied looking at the density-flux fundamental diagram 
[24], that is the behavior of the average flux as a function of the load of the network, i.e. 
the density v. The stationary distribution completely factorizes for this system, using 
the detailed balance. We consider regular lattice and homogeneous rates rj = r = 1 Vi, 
and the distribution simplifles into: 

i. V{n) =P(0)f^ if n < ni 

tt. V{n) = V{0)^^^ if ni<n<n* 

Hi. V{n) = ii n > n*. 

From the normalization condition G{1) = 1 we get 'P(O), while A can be obtained 
numerically solving with respect to A the expression for the average density G'{1) = v. 
The average flux = — V{n*)) and the average velocity v = (p/u can be easily 

computed as function of the parameters of the model. 

The behavior of the flux 0(z/) and the average velocity f (z/) are reported in Fig. 
15, for two different values of ni and for flxed capacity n* = 10. Studying these two 
quantities as a function of the particles density u, we recover several stylized facts of 
traffic flows in trasportation network [24]. The flow increases almost linearly at low 
densities and reaches a maximum when queues are half-fllled. Then the curve decreases 
until it vanishes at n*. Correspondingly, the velocity decreases, from an initial plateau 
ai V = 1 (free-motion), in a non-linear way (depending on n* and tt-i) until it vanishes at 
u = n*, where the system becomes completely overloaded. The congested phase arises 
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Figure 15. Left: Fundamental diagram for tiie model onto a regular lattice, ni = 2 
and 3, n* = 10 by analytical calculations. 

Right: Velocity profile as a function of the density for the model onto a regular lattice 
ni — 2 and 5, n* = 10 by analytical calculations. 

always continuously and no phase transition is observed. It would be interesting to 
understand if a discontinuous transition could be observed in the fundamental diagram 
by varying some of the parameters of the model. 

5. Conclusions 

We have proposed a minimal traffic model to study congestion phenomena on complex 
networks, with applications to realistic processes in information as well as transportation 
networks. 

This traffic model considers particles performing random walks through the network 
and forming queues on the nodes they visit. Particles are created on each node with 
a given rate p, but they can be absorbed (at rate fi) only during the hopping process, 
not when they are waiting into the queues. This simple rule is sufficient to generate a 
transition from a free- flow to a congested phase as a function of the creation/absorption 
rates. Moreover, the introduction of a rejection probability ri{n) for nodes with more 
than a threshold value n* of particles can modify the nature of the transition from a 
continuous to a discontinuous one. 

By means of a mean-fleld approach we have analyzed the behavior of the model on 
ensembles of generalized random graphs, where the graph properties are specifled only 
by the degree distribution P{k). The interplay between the queueing properties of the 
model and the topological structure of the network generates a rich phenomenology of 
continuous and discontinuous phase transitions. In particular, on highly heterogeneous 
networks, the system can present also an hybrid behavior in which a primary continuous 
transition to a partially congested state is followed (at higher values of the control 
parameter p) by a discontinuous transition to a state with higher level of congestion 
(possibly a completely congested state). For some explicative case we have computed 
the whole phase diagram with all possible transitions as a function of the parameters of 
the model. In general on scale free networks the critical line goes to zero in the inflnite 
size limit, but it may converge to a constant if the rejection probability is sufficiently 
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large. 

The present model is reminiscent of a class of non-equilibrium systems (e.g. Zero- 
Range Processes) in which particles condensation is observed [13]. We have shown 
that a conserved version of our model, without particles' creation and absorption, does 
present condensation and we have discussed its relation with the congestion phenomenon 
occurring in the non-conserved model. 

We have briefly described two possible apphcations of the model to packet-based 
traffic in information networks, such as the Internet at the router level, and to the 
vehicular traffic in road networks. 

In the case of information networks, the model offers a simplified picture of routing 
system, in which packets are forwarded towards destination with a routing protocol that 
can take into account the level of traffic in the local neighborhood. Our results show 
that traffic-aware routing is useful only in heterogeneous networks, where it expands the 
region of stability of the congestion-free state. However, a congested phase may arise 
abruptly, and then it may persists even under lower traffic loads (hysteresis effects). The 
model provides a simple explanation for the so-called " scahng breakdown" of the traffic 
fluctuations observed on the Internet [25], in terms of the intrinsic statistical property 
of the queues. 

On the other hand a slightly modifled version of the conserved model, in which nodes 
have a maximum capacity, can reproduce qualitatively the form of the fundamental 
diagram and velocity profile observed for the vehicular traffic in transportation networks 
[24]. 

These results show that a simple model of interacting random walks with some 
further ingredient is enough to obtain a very good description of the statistical properties 
of the collective behavior of the system. 

The work presented here can be extended in several interesting directions. The 
possibility to solve the model on a given network (e.g. Internet's map or urban road 
networks) using Eqs. 7 with realistic parameters, could provide both specific predictions 
on the robustness of the network to traffic overloads and important hints for the design 
of systems that could be less vulnerable to congestion phenomena. The dynamical 
environment created within this model could also be exploited as a framework for 
testing the statistical properties of single particle dynamics under more complex routing 
schemes, in a way similar to the study of a tracing-particle in hydrodynamics. 

Finally, it would be interesting to model the complex adaptive behavior of human 
users in communication networks, such as the Internet, by introducing variable rates of 
packets production in response to network performances. It is known that users face the 
social dilemma of maximizing their own communiction rates, maintaining the system 
far from the congested state [28]. In such a situation, the presence of a continuous 
transition may allow the system to self-organize at the edge of criticality, whereas a 
discontinuous transition may have catastrophic consequences. 
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Appendix A. Validity of Product-measure distributions 

The existence of factorized stationary states in interacting particles systems on graphs 
is a subject of great interest, because models with this property are exactly solvable 
and can be used to investigate complex dynamical phenomena like mass-transport or 
condensation. The most studied family of models with a product-measure stationary 
state are the zero-range processes (ZRP), in which it comes directly from the fact that 
hop rates depend only on the number of particles in the departure node [13]. In the 
Misanthrope processes, instead, the hop rate depends on both departure and arrival 
node, but under some conditions the stationary state distribution is still factorizable on 
single nodes [13]. More general sufficient conditions to have factorized stationary states 
in continuous mass-transport models on general graphs were pointed out recently by 
Evans et al. [29]. 

All statistical physics models on which factorization has been proved have either 
a conserved number of particles or zero-range interaction, and satisfy either a detailed 
balance condition or a pairwise balance condition, depending on whether the transition 
rates are symmetric or not along the links of the network. A partial generalization 
of these results can be obtained applying Jackson's Theorem, a well-known result in 
queucing theory, that requires only the stationarity condition and a local conservation 
law for the fluxes of particles [15]. To our knowledge there is no general result for 
factorization in models with non-conserved number of particles in which the hop rates 
depends on the number of particles in both departure and arrival nodes. The aim of 
this appendix is to discuss in detail some special cases in which factorization is exact. 

Appendix A. 1. Detailed Balance Condition 

Let us consider the conserved model in which Af particles are deployed on a general 
graph and can hop from node i to j with hop rate Uij{ni, rij), i.e. a general Misanthrope 
process. In our model, we will specify later Uij{ni, rij) — riWijll — r]j{nj)] with Wij — p. 
When there are no currents in the graph, the detailed balance should be valid, and this 
is the case of models in which the hop rate is symmetric along the edges of a node. The 
dynamics is based on random walks, thus we expect detailed balance to hold. 

Given the states n = (ni, .... ... , Uj, . . . , njv) and n/ = {ui, . . . , nj + 1, . . . , — 

1, . . . jUn), with transition rates w{n n/) = Uji{nj,ni) and w{n/ n) = Ui.j{ni + 
1, Tij — 1), imposing the detailed balance condition and introducing the factorized ansatz 
V{ni, . . . , un) = Yli 'Pii^i)^ we get the equation 



The case Ui = 0, Uj > gives Uij{l, Uj — l)T'i{l)Vj{nj — 1) = Uji{nj, 0)Vi{0)Vj{nj) , that 
is a recurrence equation producing 



Uij{ni + 1, Uj - l)Vi{ni + l)Vj{nj - 1) = Uj,{nj,ni)Vi{ni)Vj{nj) (A.l) 




(A.2) 
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For this result to be correct also for > 0, we plug it into Eq. A.l and obtain a 
condition on the hop rates 



Uij{ni + l,nj - 1) 



Uij{l,0) Uji{l,ni)uj,{nj,0) 



u 



1,0) Uij{ni + l,0)uij{l,nj - 1) 



Uji{nj,ni) (A. 3) 



The condition is satisfied by the hop rates Uij{ni, nj) — of conserved 

model, that thus admits a product-measure stationary distribution. Reasonably 

assuming r]{0) = and using properties from Eq. A. 2, the single-site distribution 
becomes 



[1 - - 1)] 



(A.4) 



where Vi{0) is fixed by the normalization and A is a constant independent of the node. 

When we introduce particle creation and destruction, detailed balance is not always 
satisfied on general graphs. An example is provided by the non-conserved model 
without rejection probability. Suppose the particles are created on node i with rate 
Pi, and absorbed when they move to j with probability /ij. Note that the particles 
transfer occurs now with rate Uij{ni,nj) — ^{1 — /ij). Writing the detailed balance 
for the particle-transfer transition and factorizing the distribution we find the relation 
|(1 - i^j)Vi{ni + l)Vj{nj - 1) = g(l - iii)Vi{ni)Vj{nj) that gives 

'rikjViiiy 



I.e. 



ki 

nA 



(A.5) 



(A.6) 



where A is again a constant independent of the node properties. This expression satisfies 
also the condition for general Ui > 0, but when plugged it into the detailed balance 
condition for the creation-absorption transition, i.e. PiVi{ni) = ^ ^j^v{i) + 1)) 

produces a node-dependent expression for A. Only in the special case of an homogeneous 
set of parameters on a regular graph (i.e. ki = k Vi), detailed balance is enough to have 
a product-measure, that reads V{n) = (1 — p/rfj.){p/rfj.)"'. 

Hence, detailed balance is not a sufficient condition for factorization of our model 
on a general graph, not even in the simpler case in which particles are not rejected 
by destination nodes. In the next section, we will check if a product-measure can be 
derived from the weaker condition imposed by stationarity. 



Appendix A. 2. Stationarity condition and the Jackson's Theorem 

The reason of the failure of detailed balance is not the particular mechanism of particles 
destruction that we have considered. Indeed, a non-conserved model with standard 
birth-death process at the nodes (i.e. particles are created with rate Pi and destroyed 
with rate /ij) would have the same problems. Nonetheless, such a model admits a 
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product-measure stationary distribution. This can be easily proved using Jackson's 
Theorem, one of the major results in queueing theory [15]. 

Jackson's Theorem is usually formulated for networks of queues with exponentially 
distributed service times and Poisson arrivals from outside and for this reason it can be 
easily adapted to describe the present system. The first ingredient of Jackson's approach 
is the global balance, or stationarity condition, of the probability distribution, 

N N N \ 

+ I] n 5] - r)j{n,)] V{m, . . . , riAr) = (A.7) 
j=i i=i j=i j 

N 



^Pi{l - Sni,o)V{ni, . . . , ni - 1, . . . , niv) 



i=l 

N N 

+ n HjWij[l - rij{nj)]V{ni, . . . , + 1, . . . , n^, . . . , njv) 
i=i j=i 

N N 

+ ~ K,o){'^ - l^i) [1 - Vii^i - i)]V{ni, . . . , - 1, . . . , rij- + 1, . . . , njv) 

i=i j=i 

with Wij = 1/ki. Wc focus here on the dynamics without rejection {rj{n) = Vn) and 
consider two important additional relations. In the stationary un-congested state, the 
incoming flux of particles in the network has to be balanced by the outcoming flux. Let 
us call Xi the stationary rate of particles entering node i, then the balance of particles 
entering and leaving the network is given by 

i i j 

The rates {Xi} are deflned by means of a local conservation law at the nodes, 

Xi=pi + {l-iii)Y,^jWji. (A.9) 

j 

If a solution exists, this system of linear equations can be solved to find {Aj} for any 
given graph and set of parameters {pi, /Xj}. 

Once we have computed {Xi}, the stationary probability distribution (in the 
uncongested phase) can be expressed in the following product-measure form, 

V{n„ ...,n^) = nn(n.) = H " 7) (^T (^-10) 

i—i i \ * / \ « / 

The proof of factorization derives straightforwardly from the stationarity condition. 
Inserting expression A.IO into Eq. A.7 we get 



13 
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i.e. 



E^ + E^^^ + E^(^'-^^o = E(^^^+^o (A.12) 



where we have used Eqs. A. 8 and A. 9. 



Appendix B. Pseudo- factorization in the general non-conserved model 

In this Appendix, we show that the probabihty distribution of the general model 
proposed in this paper is not exactly factorizable over the nodes, i.e. a product-measure 
stationary state does not exist. Nevertheless, the way in which factorization is broken 
is very peculiar and still allows one to obtain important information on the structure of 
the stationary distribution. 

Appendix B.l. Beyond Jackson's Theorem 

If we consider also rejection of particles at the nodes, Jackson's Theorem does not apply, 
though we can prove a theorem that specifies the form of the stationary distribution. 

Theorem 1 Consider an open particles system on a general graph with adjacency 
matrix {aij} in which particles enter the system at every node i with rate pi, hop from 
node i to node j with rate Uij{ni, rij) = ViWij^l — —rij{nj)] (with Wij = aij/ki), and 
leave the system with rate Uio = ri Wij/jj[l — i]j{nj)]. Let {Aj(nj,n_j)} be functions 
of the configuration {ui, = n = (^i, • • • , ^i, • • • , ?^Ar) satisfying the system of linear 
equations 

i j 
then the stationary probability distribution of particles is 

V{n„ . . . , n^) oc n = H 11 ' ^^'^^ 

i=l i=\ i=\ * 

The proof of the theorem is straightforward and consists in inserting the factorized 
form in the stationarity condition like in Jackson's Theorem, 

E^^^ + E E ^^^^ ~ ^jK)] = 

i i j 

Vijni - 1) ^^ Vijui + 1) 

l^Pi + nWij^Xj[l - (B.3) 

E^.Viiui — 1) VAn,-. + 1) , ^ 
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i.e. the r.h.s. reads 



PiTi 



^3 



Xj{nj + l,n_j)ri 



\i{ni,n_i)rj 
Pi+[l- Tjiirii - 1)](1 - iii) X ^ji^ji^j + 1> 



n 



E 



\{ni,n_ 

i 3 

where in the last hne we have used the local flux-balance condition (B.l). The 
stationarity condition thus becomes a global balance condition X^jPi — + 
1, ri_j) Ylij fJ'jWij [1 ~ Vj i'lT'j)] ) that is verified by summing over all nodes the local balance 
conditions in Eq.B.l. 

Note that the theorem holds also for more general hop rates Uij{ni,nj) and 
for creation/destruction of particles that also depend on the number of particles in 
neighboring nodes. 



Appendix B.2. Derivation of mean- field equations and role of correlations 

The mean-field approximation consists in studying the single-node behavior, neglecting 
the correlations that could exist between difi:erent nodes. The pseudo-factorized form of 
the stationary distribution suggested by Theorem 1 says that correlations only exist at 
the level of the equations defining the set of {Aj}. These self-consistent equations could 
be solved in principle for each configuration n, but this is impossible in practice and we 
have to resort to some approximation. 

Averaging over ni^rij both sides of Eq. B.l, and replacing the rejection function 
77(71) with the corresponding rejection probability x, leads to the mean-field equation 




that is valid up to the congestion transition. 

Note that the quantity \ represents the total rate of particles outcoming from node 
i. As particles hop with rate n per node, independent from the number of particles 71^, 
we expect Aj to be proportional to and to the probability that the queue in i is not 
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C(r) 



0.5 















20 



Figure Bl. Normalized correlation in the number of packets between a node and a 
node placed at distance r in a square lattice with periodic boundary conditions of size 
L = 500, p = 0.12, n = 0.2 and fj = 0.5, n* = 10. 



empty, that we call 1 — g^. With (Aj) ~ ri{l — qi), the Eq. B.6 becomes identical to the 
condition (fi) = used in Section 3.1 to derive a system of 2A'" mean-field equations for 
the node quantities {qi,Xi} that can be solved recursively on every graph. 

When all nodes are far from a congested state, factorization is almost exact, because 
the rejection term is effectively inactive. We have verified this statement computing the 
two-points correlations C(r) = {riiTij) (||i — j\\ = r) between the queue lengths of nodes 
at distance r on a regular square lattice (see Fig. Bl). For p < p^, correlations are 
completely absent. 

Appendix C. Behavior of the system for small queueing capacity 

Most of the analytical results presented in the main text, both at the ensemble level 
and on single graphs, are obtained in the limit n* — > oo. For finite n*, the calculations 
still can be solved numerically but the overall approach becomes cumbersome and much 
less transparent. We have stressed that, except for the pathological case of the average 
density z/(p) (see Sec. 3.3), the qualitative behavior of the system does not change if we 
consider finite but large n*. It is natural to ask how far we can push this approximation, 
and if the behavior for small values of n* is still qualitatively similar to that for n* — > oo. 
For n* — 1, i.e. when each node reject particles with a probabihty fj as soon as it 
contains a particle, mean-field calculations on ensembles of random graphs with a given 
degree distribution do not present further difficulties compared to the n* — > oo case. 
Fig. CI (left) reports the behavior of the congestion order parameter p[p) on a scale-free 
network with rejection probability 77 (n) = f]6{n — 1) for f/ = 0.1 (black circles) and 
fj = 0.9 (red squares). Fig. CI (left) shows that increasing fj a discontinuous transition 
to the congested phase appears, but without any shift of the transition point to higher 
values of p. Hence, in this case, the traffic-aware protocol is not effective in enlarging 
the free phase. Nevertheless, for n* — 2 (right panel in Fig. CI) the behavior is already 
qualitatively the same that in the limit n* — > 00. 
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Figure CI. Behavior of the congestion parameter p{p) on scale- free network {N = 
3000, 7 = 3) for ^ = 0.1 (black circles), and 0.9 (red squares) with n* = 1 (left) and 
n* = 2 (right). 




Figure Dl. Time series of the total number of packets from simulations onto a square 
lattice 100X100 with pbc, ^ = 0.2, rj = 0.9, n* = 10 starting from the free phase, for 
different values of p. At p = 0.15, after a transient, the system jumps into the congested 
state 



Appendix D. Metastability of the free phase 

Since for strong enough routing protocols (high r]), in a certain range of p there is 
coexistence of the two phases, we would expect, because of finite size effects, that 
fluctuations can trigger a jump from the free to the congested phase. As we can see 
in fig. Dl this is the case for a square lattice of 100X100 nodes, in which at p = 0.15 
7] = 0.9 the system, starting from the free phase, becomes congested after a while. 
Moreover, once the system is congested into the coexistence region, the jump back into 
the free phase is not possible anymore, because the congested phase is absorbing in the 
same way of the growing phase of the pinning transition: the more the system remains 
in the congested phase, the more the queues are growing and it is more hard to turn 
back into the free phase. However, as the jump is triggered by independent activation 
processes, the transient time goes exponentially with the systems' size. 
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